Heat Conduction and Entropy Production in a One-Dimensional Hard-Particle Gas 
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We present large scale simulations for a one-dimensional chain of hard-point particles with alter- 
nating masses. We correct several claims in the recent literature based on much smaller simulations. 
Both for boundary conditions with two heat baths at different temperatures at both ends and from 
heat current autocorrelations in equilibrium we find heat conductivities k to diverge with the num- 
ber A'^ of particles. These depended very strongly on the mass ratios, and extrapolation to A'^ — » oo 
resp. f — > cxD is difficult due to very large finite-size and finite-time corrections. Nevertheless, our 
data seem compatible with a universal power law k ~ A'^" with a ~ 0.33. This suggests a relation 
to the Kardar-Parisi-Zhang model. We finally show that the hard-point gas with periodic boundary 
conditions is not chaotic in the usual sense and discuss why the system, when kept out of equilibrium, 
leads nevertheless to energy dissipation and entropy production. 



Low-dimensional systems are special in many ways. 
Second order phase transitions have anomalous expo- 
nents, chemical reactions do not follow the mass action 
law hydrodynamics breaks down due to divergent 
transport coefficients caused by long time tails ||^, and 
electrons in disordered systems are localized [|. A last 
item in this list is the divergence of heat conductivity Q 
in < 2 spatial dimensions. 

For ordered (periodic) harmonic systems it is well 
known that all transport coefficients are infinite due to 
the ballistic propagation of modes. Thus one needs ei- 
ther disorder or nonlinear effects in order to have finite 
conductivity k. For electric conductivity, disorder in 1-d 
leads to zero conduction. The main difference between 
heat and (electronic) charge conduction is that there is no 
background lattice in the former, i.e. translation invari- 
ance is not broken even if the system is disordered. Of 
course one can study the electronic contribution to heat 
conduction, but experimentally one never can neglect the 
ionic contribution [^. Thus one has always soft (Gold- 
stone) modes in heat conduction. These soft modes are 
not localized by disorder ^ , and they are not affected by 
nonlinearities. Thus they propagate essentially freely. In 
high dimensions this has no dramatic consequences. But 
in low dimensions they become important and lead to the 
above mentioned divergence. More precisely, one expects 
a power divergence k ~ A'^" in d = 1 and a logarithmic 
divergence in d = 2. Simulations and calculations using 
the Green-Kubo formula give a « 0.35 to 0.45 Q]. It is 
not clear whether the slight discrepancies found between 
different models are true or artifacts. Theoretically, one 
would of course prefer a universal value. 

There are some exceptions to this general scenario. 
Apart from models with external potentials and broken 
translation invariance, the best known 1-d model with 
finite K is the rotor model of (tJ ||] . Here very fast rotors 
effectively decouple from their less fast neighbours. Thus 
very steep jumps of temperature are effectively stabilized 
and act as barriers to energy transport. 

Another model which was recently claimed to have k fi- 
nite p is the 1-d hard point gas with alternating masses. 
The same gas with all particles having the same mass 



is trivial (a collision between particles is indistinguish- 
able from the particles just passing through each other 
undisturbed), and perturbations just propagate ballis- 
tically, leading thus to infinite k (i.e., no temperature 
gradients can build up inside the gas). To break this in- 
tegrability, it is sufficient to use alternating masses: every 
even-numbered particle has mass mi , and every odd has 
7712 = rmi with r > 1 I, IToHn], |l|, p, |l|, |l|. 

The arguments given above for a divergence of k with 
A^ hold also for the hard point gas. It is obviously nonlin- 
ear, sound waves dissipate, it is translation invariant, and 
there is no obvious special feature as in the rotor model. 
Indeed, prior to heat conduction has been studied 
in this system by means of simulations in |l^, |l^. In 
these papers it was claimed that k diverges. But the 
simulations of | p3[ had presumably low statistics (very 
few details were given), while the simulations of are 
obviously not yet in the scaling regime {N is too small) 
and are compatible with k const for N oo. In any 
case, the exponent a suggested by the simulations of |l^ 
is < 0.22, much smaller than for all other models. 

In view of this confusing situation, and suspecting that 
the simulations of |l^ Q were not done most effi- 
ciently, we decided to make some longer simulations. 

We followed ||l^ in setting ttii — 1 and using 
Maxwellian heat baths at the ends with Ti = 2,T2 = 8 
(i.e., after hitting the end, a particle is reflected with 
a random velocity distributed according to P{v) = 
e{±v)m,v/TexY>{-mv^/2Ti^2)- The heat baths sit at 
a; = and x — N, i.e. the density of the gas is 1. When 
an even particle with velocity vi collides with an odd one 
with velocity V2, their velocities after the collision are 



(1 — r)vi + 2rv2 
(1 + 



27;i + (r - l)v2 
(1 + r) 



(1) 



Between two collisions, the particles propagate freely. 
Thus a fast simulation algorithm is event-driven: For 
each particle i we remember its velocity, the time ti of 
its last collision (initially, ti is put to zero), and the po- 
sition Xi it had at that time. In addition, we maintain a 
list of future collision times Ti for each neighboring pair 
+ 1) (here the walls are treated formally as parti- 
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cles with vo = vn+i = 0, Xq = 0, Xn+i = N). The 
system is evolved by searching the smahest r^, updating 
the triples {ti,Xi,Vi) and (i^+i, x^+i, Wi+i), and calculat- 
ing the new future collision times Ti_i and r^+i (the new 
Ti is infinite). Since the list {ri} is essentially a priority 
queue jl^, we use for it the appropriate data structure 
of a heap Using heaps, searching for the next colli- 
sion takes a CPU time O(lnA^). In comparison, a naive 
search would take 0{N). This allowed us to make much 
larger simulations than in previous works. Our largest 
systems contained 16383 particles and were followed for 
> 10^^ collisions. In spite of this, we had to start with 
carefully tailored initial configurations to keep transients 
short. When obtaining statistics one should not forget 
that measurements should not be made after a fixed num- 
ber of collisions, but at fixed intervals in real time. The 
correctness of the results and the absence of transients 
were checked by verifying that the energy density is con- 
stant, as proven in ||l4| . 

Before presenting results, let us discuss the expected 
dependence on the mass ratio r. For r — s- 1, equilibration 
becomes slow (it takes a long time until a fast particle is 
slowed down to average speed), but perturbations prop- 
agate ballistically. Thus a perturbation will be damped 
out slowly at first, but it will have no long time tails 
and is damped exponentially. In the other extreme, for 
r — > oo the light particles bounce between pairs of heavy 
ones which are hardly perturbed. Thus, if a heavy parti- 
cle is perturbed, we have a situation very similar to the 
one for r — > 1. If a light particle is perturbed, its en- 
ergy is soon given to its two nearest heavy neighbours, 
which then behave again as for r « 1. In contrast, in 
the intermediate region 1 <C r <^ cx) we expect the per- 
turbation to spread non-ballistically for all times. It is 
in this regime that we expect the fastest convergence to 
asymptotic behaviour, both with respect to time and to 
N. 

In Fig. la we show k, defined as total energy flux J = 
"^iiTiivf /2 divided by AT, versus N, for four values of 
r. The value r — 1.22 is in the small-r region and was 
studied most intensively in Q. The value r — 2.62 is 
near the center of the intermediate regime, while r — 5 
is clearly above it. Finally, r = 1.618 = (1 -I- a/5)/2 
was chosen because it was used in [|| , not because of its 
irrationality (problems with ergodicity related to rational 
values of r exist only for very small N [|ll|, |l^). 

A power law would give a straight line with slope a 
in Fig. la. None of the four curves is really straight. For 
small the curve for r — 1.22 agrees perfectly with Fig. 3 
of (which extends only to = 1281). It shows the 
strongest curvature (in agreement with the above dis- 
cussion), and the small- A^ data alone would suggest a 
cross-over to k = const. But this curvature stops for 
large A^ and a closer look shows that the slope increases 
for A^ > 8000. The same is true also for the other curves: 
They all bend down for small N but veer up for larger 
systems (Fig. lb). This is most clearly seen for r = 1.618 
and 2.62. It is less clear for r = 5, but the most ratio- 
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FIG. 1: (a); Log-log plot of J/{T2 - Ti) versus for four 
values of r. Statistical errors are always smaller than the 
data symbols, (b): Part of the same data divided by A'^" with 
a = 0.32, so the y-axis is much expanded. 

nal expectation is that also this curve will have the same 
slope for A oo. Our best estimate a — 0.32 ^g ^ij has 
asymmetric errors because we do not know how much 
more the curves will bend upwards for very large A^. 

The rescaled temperature profiles for r — 1.22 are 
shown in Fig. 2. To verify the claim of Q that T{x) 
approaches the profile Tk{x) predicted by kinetic theory 
with an inverse power of A^, i.e. T{x) — Tk{x) ^ A^^°-^^, 
we plot T{x) - Tk{x) against x/N. For A^ < 2000 we see 
indeed this convergence in perfect agreement with [ p^ , 
but not for A^ > 2000. Instead, it seems that T{x)-Tk{x) 
remains different from zero for N oo. The analogous 
results for r = 1.618 are shown in Fig. 3. In that case, 
the scaling observed in |l^ is confined to very small A^, 
not shown in the figure. The fact that T{x) — Tk{x) re- 
mains finite for A?^ ^ oo is now obvious. In contrast to 
a conjecture in ||], the temperature profile also does not 
become linear for large A^. All results for r — 1.618 are 
qualitatively also true for r = 2.62 (not shown). 

These results are easily understood. For r = 1.22 we 
are in the small-r regime. This explains the slow conver- 
gence of a with N and the weakness of long time tails, 
manifested in the agreement with kinetic theory. Only 
at very large A^ we do see the correct asymptotics. For 
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FIG. 2: T(x) - rfc(a;) against x/N for r = 1.22, where 
Tk{x) = [T^^/^ - {T'I''-'' -T'^'^)x/Nf^^ is the temperature pro- 
file according to kinetic theory. In order to reduce statistical 
fluctuations, we averaged in the curves for A'' > 1023 over 
three successive values of x. 
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FIG. 3: T{x) - Tk(x) against x/N for r = 1.618. 



r — 1.618 and 2.62 we are no longer in this regime, the 
long time tails are stronger, and the disagreement with 
kinetic theory is more obvious. 

In addition to systems driven by thermostats at dif- 
ferent temperatures, we also studied systems in equilib- 
rium with periodic boundary conditions. Here the Green- 
Kubo formula allows k to be calculated from an integral 
over the heat current autocorrelation (J(t)J(O)) Q]. In 
Fig. 4 we show this, after suitable normalization and after 
multiplication by a power of t which makes it constant 
for large N and t. We see strong oscillations with peri- 
ods oc N which reflect the dominance of (damped) sound 
waves with a fixed velocity of sound (see also |15t ) . They 
were mistaken for statistical fluctuations in |9|], showing 
clearly that the simulations of have not reached the 
asymptotic regime in contrast to what the authors as- 
sumed. Our data suggest that (J(t)J(O)) ~ t'^-^'^ for 
large N with a cut-off at t oc A^, giving a — 0.34 in 
perfect agreement with our previous estimate. 

A 1-d hard particle gas should be described macro- 
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FIG. 4: Total heat current autocorrelation, t°-'^'^N~'^ 
{J{t)J(0)) for r = 2.2 and T = 2. Total momentum is P = 0. 



scopically by hydrodynamics, i.e. by the Burgers resp. 
Kardar-Parisi-Zhang (KPZ) equation ^7^. If we assume 
that heat diffusion scales like diffusion in KPZ, we might 
expect a = 1/3 in agreement with our numerics. But 
we should warn that particle spreading in the 1-d hard 
particle gas is not superdiffusive (unpublished data; for 
r = 1 see |l^), so the relation with KPZ is not trivial. 

Finally, we measured also the propagation of infinitesi- 
mal perturbations. Similar simulations were also made in 
1^ , but there the perturbations were added to the ground 
state {E = 0, all particles are at rest). In contrast, we 
perturbed equilibrium states, i.e. we performed a stan- 
dard stability analysis as used e.g. to estimate Lyapunov 
exponents. Indeed, we only followed the perturbation in 
velocity space, not in real space. More precisely, after 
having run the system long enough to have eliminated 
transients, we chose a tangent vector 



(<5z;,(0),fe,(0)) = (5,,o,0) 



(2) 



and iterated, together with the system itself, its lin- 
earized variational equations. Notice that 6vi(t) are inde- 
pendent of the positions perturbations Sxi (t) for nearly 
all times (whenever there is no collision), thus it is pos- 
sible and legitimate to solve the variational equations for 
velocities only. After the evolution of the Svi {t) has been 
followed for a given time tmax, Svi{0) is again chosen as 
Si^o, and the integration is continued. 

It is easy to prove that that the 1-d hard point 
gas is not chaotic in the usual sense. Consider the 
weighted L2 norm of the perturbation, ||(5w(t)||2 = 
['^i^imi{Svi{t))^]^/^ . Since the Svi{t) change during 
a scattering according to the same Eq.(P as the ve- 
locities Vi(t) themselves, energy conservation leads to 
||(5w(i)||2 = 1. Indeed the absence of chaos is quite ob- 
vious since there is no local instability. It seems to con- 
tradict a widespread folklore that dissipation and entropy 
production are tightly related to chaos (which sometimes 
is true; e.g. ||l^, page 231) - although it is also appreci- 
ated that this might not be always the case pO| . 
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FIG. 5: Effective diffusion coefficients for various values of r. 
Statistical errors are smaller than the sizes of the symbols. 
Temperature is T = 2, total momentum is P = 0. 

One solution to this puzzle is the observation go- 
ing back to work by Wolfram on cellular automata [ p2| , 
that the notion of chaos in systems with infinitely many 
degrees of freedom is ambiguous and is not necessarily 
related to any local instability. In a spatially extended 
system it makes perfect sense to use a norm which, in 
contrast to the L2 norm used above, puts most weight 
on near-by regions and exponentially little weight on far- 
away regions. With such a definition, the norm of a per- 
turbation moving towards (away from) the observer with 
constant velocity will increase (decrease) exponentially. 
More generally, also perturbations spreading diffusively 
will lead to an increase of the uncertainty about the lo- 
cal state for a short-sighted observer. For the 1-d hard 
point gas this means that there is no need for any local 
instability to generate dissipation, local thermal equilib- 
rium, and mixing. In a non-equilibrium case entropy flow 
is provided by the stochastic thermostats at the ends. 



while (coarse-grained) entropy is produced by the diffu- 
sive propagation of perturbations. 

Indeed, the situation is not quite as simple due to the 
divergence of k with N which suggests that perturbations 
propagate super-diffusively. For a more direct proof, 
we show in Fig. 5 their effective diffusion coefficient, i.e. 
the average squared distances divided by t (notice that 
^m,(<5«.(t))2 = l), 

N 

D = {X^)/t = i^rn,{5vdt)f . (3) 

For diffusive spreading, this would be constant. Instead, 
it increases with t for all r. This increase does not follow 
a pure power law, but again deviations from a power law 
are strongest for very small and very large r. For mod- 
erately large t they would suggest a crossover to normal 
diffusion {D — const), but for all r > 1.6 this is ruled out 
by the values at very large t. Again it is hard to estimate 
the asymptotic behaviour precisely, because the curves 
bend upward for large t. If we assume we obtain 

a precise lower bound (3 > 0.36 ±0.01, but only a poor 
upper bound which just excludes ballistic behaviour. 

In summary, we have given compelling evidence that 
heat conduction in the 1-d hard point gas shows the 
anomalous divergence with system size expected for any 
generic 1-d system, in spite of strong finite-size and finite- 
time effects. This "normal" anomalous behaviour holds 
in spite of the fact that the system is not chaotic in the 
usual sense, proving again that chaos in the form of local 
instabilities is not needed for mixing behaviour and dissi- 
pation. Finally, we have discussed a possible connection 
to KPZ scaling. 

We are indebted to Roberto Livi and Antonio Politi 
for very helpful correspondence. W. N. is supported by 
the DFG, Sonderforschungsbereich 237. 
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